program BuoyCollocate
!
! Program to collocate buoy data with scatterometer winds.
! The buoy data must be in ASCII format as written by BuoyRead
! Output data are also written in ASCII format
!
! Last Modified on: $Date: 2016-05-31 11:23:13 +0200 (Tue, 31 May 2016) $
!

  !  #[ modules used
  use BufrMod, only: set_BUFR_fileattributes, get_BUFR_nr_of_messages, &
                     open_BUFR_file, get_BUFR_message, close_BUFR_file, &
                     bufr_file_attr_data, BufrDataType, &
                     missing_real, missing_int
  use BufrMod, only: missing_indicator_real_r4=>missing_indicator_real

  use bufr_data_extraction_module, only:  set_bufr_template, goto_subset, &
       get_lat, get_lon, get_qcflag, get_res, get_sw_id_l1, &
       get_year, get_month, get_day, get_hour, get_minute, &
       get_dir, get_spd, get_mdir, get_mspd, get_cln, get_trk, get_bsdist, &
       get_aspd, get_adir

  use Compiler_Features, only: getarg_genscat, iargc_genscat
  use DateTimeMod, only: ymd2julian
  use convert, only: get_distance
  use ErrorHandler, only: no_error
  !  #]
  !  #[ variables and types
  implicit none

  ! data type for collocations
  type :: buoy_coll_type
    ! data read from the buoy ascii input file
    integer :: buoy_id
    real    :: buoy_lat
    real    :: buoy_lon
    integer :: buoy_year
    integer :: buoy_month
    integer :: buoy_day
    integer :: buoy_hour
    integer :: buoy_minute
    real    :: buoy_speed
    real    :: buoy_dir
    real    :: buoy_temp
    real    :: buoy_prs
    real    :: buoy_rhu
    ! this one stores the julian date, generated from the above items
    real    :: buoy_julian
    ! data read from the scat bufr input file
    real    :: scat_lat
    real    :: scat_lon
    integer :: scat_year
    integer :: scat_month
    integer :: scat_day
    integer :: scat_hour
    integer :: scat_minute
    real    :: scat_speed
    real    :: scat_dir
    real    :: model_speed
    real    :: model_dir
    real    :: distance
    integer :: wvc
    real    :: alfa
    real    :: mle
    logical :: qcflag
    ! optional inputs read from the scat bufr input file
    ! only used when the research_option_set is .true.
    real    :: analysis_speed 
    real    :: analysis_dir
  end type buoy_coll_type
    
  type(buoy_coll_type), dimension(:), allocatable :: buoy_coll

  type(bufr_file_attr_data) :: bufr_file_attributes
  type(BufrDataType)        :: tbd
  integer            :: kmsg, ksub, nr_of_messages
  integer            :: kbuoy, nbuoys
  integer            :: wvc, res
  integer            :: narg, iarg, error_flag
  integer            :: year, month, day, hour, minute
  integer            :: jul_tmp, julian_reference
  real               :: lat, lon
  real               :: spd, dir, jul, mspd, mdir, alfa, mle, aspd, adir
  real               :: buoy_speed, buoy_prs, buoy_rhu, buoy_temp, ps, pv, pd, rho
  real               :: distance, max_distance, max_dist_deg
  real               :: time = 30.0
  logical            :: qcflag
  logical            :: print_rejected = .false.
  logical            :: print_model    = .false.
  logical            :: stresseq       = .false.
  character(len=256) :: filename_bufr, filename_output
  character(len=256) :: filename_buoy, filename_ascii
  character(len=256) :: arg, BUFR_TABLES, tmpstr
  character(len=20)  :: dtg1, dtg2
  integer, save      :: stored_res

  integer,parameter  :: unit_ascii_output     = 29
  integer,parameter  :: unit_buoy_ascii_input = 30
  integer,parameter  :: unit_scat_ascii_input = 31
  logical            :: research_option_set = .false.
  logical            :: end_of_file_reached
  logical            :: use_ascii_scat_input, use_bufr_scat_input

  ! in this program, we convert all data from the BUFR module to default
  ! size reals (4 bytes, generally)
  ! to prevent rounding errors when checking for missing values, the
  ! value of missing_indicator_real is also rounded
  real, parameter :: missing_indicator_real = real(missing_indicator_real_r4)

  ! needed for interpretation of qcflags
  integer, parameter :: instr_undefined   = -1
  integer, parameter :: instr_is_ascat    = 1
  integer, parameter :: instr_is_seawinds = 2
  integer, save :: current_instrument = instr_undefined
  integer :: instrument
  !  #]
  !  #[ initialisations
  filename_bufr(:) = ' '
  filename_ascii(:) = ' '
  filename_output(:) = ' '
  filename_buoy(:) = ' '
  res    = -1
  wvc    = -1
  lat    = 999.
  lon    = 999.
  mspd   = 999.
  mdir   = 999.
  qcflag = .false.
  alfa   = 999.
  aspd   = 999.
  adir   = 999.
  spd    = 999.
  dir    = 999.
  mle    = 999.
  year   = -1
  month  = -1
  day    = -1
  hour   = -1
  minute = -1
  research_option_set = .false.
  !  #]
  !  #[ commandline handling

  call GetEnv('BUFR_TABLES', BUFR_TABLES)   ! Environment variable
  if (len_trim(BUFR_TABLES) .eq. 0) then
    print '(A)','Please set ${BUFR_TABLES}! Processing halted!'
    stop
  endif

  ! read command line
  narg = iargc_genscat()
  if (narg .eq. 0) then
    print '(A)',' '
    print '(A)','Usage: BuoyCollocate -b <buoy ascii input data file>'
    print '(A)','                     -o <ASCII output file with colloc result>'
    print '(A)','       [and only one of the next two inputs]'
    print '(A)','                     -f <BUFR input file with scat data>'
    print '(A)','                     -a <ASCII input file with scat data>'
    print '(A)','                    [-t <max time interval in minutes>]'
    print '(A)','                    [-rejected] [-model] [-stresseq]'
    print '(A)',' '
    stop
  endif

  iarg = 0
  do while (iarg .lt. narg)
    iarg = iarg + 1
    call getarg_genscat(iarg,arg)

    select case(arg)

    ! bufr input file with scat data
    case('-f')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_bufr)

   case('-a')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_ascii)

    ! colloc result ascii output file
    case('-o')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_output)

    ! buoy ascii input data file
    case('-b')
      iarg = iarg + 1
      call getarg_genscat(iarg,filename_buoy)

    ! time interval in minutes between collocations
    case('-t')
      iarg = iarg + 1
      call getarg_genscat(iarg,tmpstr)
      read(tmpstr,fmt=*) time

    ! scatterometer or model data
    case('-model')
      print_model = .true.

    ! accepted or rejected data
    case('-rejected')
      print_rejected = .true.

    ! accepted or rejected data
    case('-stresseq')
      stresseq = .true.
    end select
  enddo

  ! sanity check
  if (trim(filename_bufr) .eq. '') then
     if (trim(filename_ascii) .eq. '') then
        print *,"ERROR: one of the 2 options -a or -f must be used."
        stop
     end if
  end if
  if (trim(filename_bufr) .ne. '') then
     if (trim(filename_ascii) .ne. '') then
        print *,"ERROR: only one of the 2 options -a or -f may be used."
        stop
     end if
  end if

  use_bufr_scat_input = .false.
  use_ascii_scat_input = .false.

  if (trim(filename_bufr) .ne. '') use_bufr_scat_input=.true.
  if (trim(filename_ascii) .ne. '') use_ascii_scat_input=.true.

  !  #]
  !  #[ read buoy data
  ! count number of buoy data in file
  error_flag = 0
  nbuoys = 0
  open(unit_buoy_ascii_input, file=trim(filename_buoy), status='old')
  do while (error_flag .eq. 0)
    read(unit_buoy_ascii_input, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0 .and. tmpstr(1:1) .ne. '#') then
      nbuoys = nbuoys + 1
    endif
  enddo
  close(unit_buoy_ascii_input)

  allocate(buoy_coll(nbuoys), stat=error_flag)
  if (error_flag .ne. 0) then
    print '(A)',' '
    print '(A,I4)','Error in allocation of buoy_coll, status =', error_flag
    print '(A)',' '
    stop
  endif

  ! read buoy data in buoy_coll structure
  ! each buoy measurement can yield 0 or 1 collocations
  nbuoys = 0
  julian_reference = -1
  open(unit_buoy_ascii_input, file=trim(filename_buoy), status='old')
  do while (error_flag .eq. 0 .and. nbuoys .lt. size(buoy_coll))
    read(unit_buoy_ascii_input, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0 .and. tmpstr(1:1) .ne. '#') then

      ! read formatted file, stop if format seems to be invalid
      nbuoys = nbuoys + 1
      read(tmpstr, fmt='(I5,F7.2,F8.2,12X,I4,2I2,1X,2I2,F7.1,7X,F8.0,F4.0,6X,F7.1,F6.1)', iostat=error_flag) &
        buoy_coll(nbuoys)%buoy_id, &
        buoy_coll(nbuoys)%buoy_lat, buoy_coll(nbuoys)%buoy_lon, &
        buoy_coll(nbuoys)%buoy_year, buoy_coll(nbuoys)%buoy_month, buoy_coll(nbuoys)%buoy_day, &
        buoy_coll(nbuoys)%buoy_hour, buoy_coll(nbuoys)%buoy_minute, &
        buoy_coll(nbuoys)%buoy_temp, buoy_coll(nbuoys)%buoy_prs, buoy_coll(nbuoys)%buoy_rhu, &
        buoy_coll(nbuoys)%buoy_dir, buoy_coll(nbuoys)%buoy_speed

      if (error_flag .ne. 0) then
        print '(A)',' '
        print '(3A,I6)','File ', trim(filename_buoy), &
          ' seems to contain invalid buoy data near line', nbuoys
        print '(A)',' '
        stop
      endif

      ! calculate julian day
      call ymd2julian(buoy_coll(nbuoys)%buoy_year, buoy_coll(nbuoys)%buoy_month, buoy_coll(nbuoys)%buoy_day, &
                      jul_tmp)

      ! subtract a reference julian day in order to prevent rounding problems
      ! since integers for julian dates are too large to be properly
      ! represented by 4 byte reals.
      if (julian_reference .eq. -1) then
        julian_reference = jul_tmp
      endif
      jul_tmp = jul_tmp - julian_reference

      ! store the julian day and the day fraction (hours/minutes)
      buoy_coll(nbuoys)%buoy_julian = real(jul_tmp) + real(buoy_coll(nbuoys)%buoy_hour) / 24.0 + &
                                                      real(buoy_coll(nbuoys)%buoy_minute) / 1440.0
      buoy_coll(nbuoys)%distance = missing_indicator_real
    endif
  enddo
  close(unit_buoy_ascii_input)
  !  #]

  if (use_bufr_scat_input) then
     call open_scat_bufr_file(filename_bufr)
  end if

  if (use_ascii_scat_input) then
     instrument = instr_is_ascat ! hardcoded for now
     call open_scat_ascii_file(filename_ascii, instrument)
  end if

  ! loop over messages in BUFR file
  max_distance = missing_indicator_real
  max_dist_deg = missing_indicator_real

  end_of_file_reached = .false.
  scatdataloop: do while (.not. end_of_file_reached)

     if (use_bufr_scat_input) then
        call get_next_bufr_scat_result(res, wvc, lat, lon, &
             mspd, mdir, qcflag, alfa, aspd, adir, spd, dir, &
             mle, year, month, day, hour, minute, &
             research_option_set, end_of_file_reached)
     end if

     if (use_ascii_scat_input) then
        call get_next_ascii_scat_result(res, wvc, lat, lon, &
             mspd, mdir, qcflag, alfa, aspd, adir, spd, dir, &
             mle, year, month, day, hour, minute, &
             research_option_set, end_of_file_reached)
     end if

     if (end_of_file_reached) exit scatdataloop
     
     if (missing_real(lat)  .or. missing_real(lon))  cycle scatdataloop
     if (missing_real(mspd) .or. missing_real(mdir)) cycle scatdataloop
     if (missing_real(spd)  .or. missing_real(dir))  cycle scatdataloop
     if (missing_real(mle)) cycle scatdataloop

     if (missing_real(max_distance)) then
        if (.not. missing_int(res)) then
           ! calculate max distance between buoy and WVC in km 
           ! and degrees (assume 1 km = 0.01 degree)
           ! max distance is cell spacing devided by sqrt(2)
           max_distance = real(res) * 0.7071 / 1000.0
           max_dist_deg = max_distance * 0.01
        endif
     endif

     if (print_rejected) then
       if (.not. qcflag) cycle scatdataloop  ! print rejected data
     else
       if (qcflag) cycle scatdataloop        ! print accepted data
     endif

     ! calculate julian date and subtract the same offset as in the buoy winds
     if (.not. missing_int(year)) then
        call ymd2julian(year, month, day, jul_tmp)
        jul_tmp = jul_tmp - julian_reference
        jul = real(jul_tmp) + &
              real(hour) / 24.0 + &
              real(minute) / 1440.0
     else
        cycle scatdataloop
     endif

     ! calculate satellite direction of flight
     ! counterclockwise with 0 degrees pointing to the North
     alfa = 360.0 - alfa
     if (alfa .lt. 0.0 .or. alfa .gt. 360.0) then
       alfa = 999.0
     endif
     
     ! step through list of buoy data
     collloop:do kbuoy = 1, nbuoys

!if (abs(jul - buoy_coll(kbuoy)%buoy_julian) .lt. 30.0/1440.0) then
!   if (abs(lat - buoy_coll(kbuoy)%buoy_lat) .lt. max_dist_deg) then
!      print *,'DEBUG: jul1',jul 
!      print *,'DEBUG: jul2',buoy_coll(kbuoy)%buoy_julian
!      print *,'DEBUG: lat', lat
!      print *,'buoy_coll(kbuoy)%buoy_lat = ', buoy_coll(kbuoy)%buoy_lat
!   end if
!end if


        ! first collocation selection criteria: time difference < specified max and
        ! lat difference < max_dist_deg (reduce nr of calls to get_distance)
        if (abs(jul - buoy_coll(kbuoy)%buoy_julian) .lt. time/1440.0 .and. &
            abs(lat - buoy_coll(kbuoy)%buoy_lat) .lt. max_dist_deg) then

          ! calculate distance between buoy and WVC
          distance = get_distance(lat, lon, buoy_coll(kbuoy)%buoy_lat,&
                                  buoy_coll(kbuoy)%buoy_lon)

          ! we have a collocation if the distance is small enough
          ! if we already had one and the current WVC is closer,
          ! the collocation is replaced
          if ( (distance .lt. max_distance)             .and. &
               (distance .lt. buoy_coll(kbuoy)%distance)       ) then
            buoy_coll(kbuoy)%distance    = distance
            buoy_coll(kbuoy)%scat_lat    = lat
            buoy_coll(kbuoy)%scat_lon    = lon
            buoy_coll(kbuoy)%scat_year   = year
            buoy_coll(kbuoy)%scat_month  = month
            buoy_coll(kbuoy)%scat_day    = day
            buoy_coll(kbuoy)%scat_hour   = hour
            buoy_coll(kbuoy)%scat_minute = minute
            if (print_model) then
              buoy_coll(kbuoy)%scat_speed  = mspd  ! print model data
              buoy_coll(kbuoy)%scat_dir    = mdir  ! print model data
            else
              buoy_coll(kbuoy)%scat_speed  = spd
              buoy_coll(kbuoy)%scat_dir    = dir
            endif
            buoy_coll(kbuoy)%model_speed = mspd
            buoy_coll(kbuoy)%model_dir   = mdir
            buoy_coll(kbuoy)%wvc         = wvc
            buoy_coll(kbuoy)%alfa        = alfa
            buoy_coll(kbuoy)%mle         = mle
            buoy_coll(kbuoy)%qcflag      = qcflag
            if (research_option_set) then
              buoy_coll(kbuoy)%analysis_speed = aspd
              buoy_coll(kbuoy)%analysis_dir   = adir
            end if
          endif
        endif
     enddo collloop

   enddo scatdataloop

  !  #[ write the ascii output file

  ! open ASCII output file and write heading, unless no collocations were found
  open(unit_ascii_output, file=filename_output)

  if (research_option_set) then
    write(unit_ascii_output,'(A)') &
    '#Buoy    Lat     Lon      Date Time Speed   Dir Scat Lat     Lon' //  &
    '      Date Time Speed   Dir  MSpd  MDir WVC    Alfa     MLE  Aspd  Adir'
  else
    write(unit_ascii_output,'(A)') &
    '#Buoy    Lat     Lon      Date Time Speed   Dir Scat Lat     Lon' //  &
    '      Date Time Speed   Dir  MSpd  MDir WVC    Alfa     MLE'
  end if

  ! write the collocations
  do kbuoy = 1, nbuoys
    if (.not. missing_real(buoy_coll(kbuoy)%distance)) then
      write(dtg1,101) buoy_coll(kbuoy)%buoy_year, buoy_coll(kbuoy)%buoy_month, &
                      buoy_coll(kbuoy)%buoy_day, &
                      buoy_coll(kbuoy)%buoy_hour, buoy_coll(kbuoy)%buoy_minute
      write(dtg2,101) buoy_coll(kbuoy)%scat_year, buoy_coll(kbuoy)%scat_month, &
                      buoy_coll(kbuoy)%scat_day, &
                      buoy_coll(kbuoy)%scat_hour, buoy_coll(kbuoy)%scat_minute

      ! convert equivalent neutral wind to stress equivalent wind if requested
      buoy_speed = buoy_coll(kbuoy)%buoy_speed
      if (stresseq) then
        buoy_temp = buoy_coll(kbuoy)%buoy_temp
        if (buoy_temp .gt. 233.16 .and. buoy_temp .lt. 313.16) then

          ! valid temperature, convert to stress equivalent wind
          buoy_prs = buoy_coll(kbuoy)%buoy_prs
          if (buoy_prs .le. 84000.0 .or. buoy_prs .ge. 105000.0) then

            ! no valid pressure, take average pressure value over 2015
            buoy_prs = 101467.0
          endif
          buoy_rhu = buoy_coll(kbuoy)%buoy_rhu * 0.01  ! percentage to fraction
          if (buoy_rhu .le. 0.0 .or. buoy_rhu .ge. 1.0) then

            ! no valid relative humidity, take average rhu over 2015
            buoy_rhu = 0.80
          endif

          ! compute saturation vapour pressure according to ECMWF IFS method
          ps = 611.21 * exp(17.502 * (buoy_temp - 273.16) / (buoy_temp - 32.19))

          pv = buoy_rhu * ps  ! actual water vapour pressure
          pd = buoy_prs - pv  ! actual dry air pressure

          ! air density is sum of partial pressures times their molar weights
          ! divided by the universal gas constant times the temperature
          rho = (pd * 0.028964 + pv * 0.018016) / (8.314 * buoy_temp)

          buoy_speed = buoy_coll(kbuoy)%buoy_speed * sqrt(rho / 1.225)
        endif
      endif

      if (research_option_set) then
        write(unit_ascii_output,103) buoy_coll(kbuoy)%buoy_id, &
          buoy_coll(kbuoy)%buoy_lat, buoy_coll(kbuoy)%buoy_lon, &
          trim(dtg1), buoy_speed, buoy_coll(kbuoy)%buoy_dir, &
          buoy_coll(kbuoy)%scat_lat, buoy_coll(kbuoy)%scat_lon, trim(dtg2), &
          buoy_coll(kbuoy)%scat_speed, buoy_coll(kbuoy)%scat_dir, &
          buoy_coll(kbuoy)%model_speed, buoy_coll(kbuoy)%model_dir, &
          buoy_coll(kbuoy)%wvc, buoy_coll(kbuoy)%alfa,  buoy_coll(kbuoy)%mle, &
          buoy_coll(kbuoy)%analysis_speed, buoy_coll(kbuoy)%analysis_dir
      else
        write(unit_ascii_output,102) buoy_coll(kbuoy)%buoy_id, &
          buoy_coll(kbuoy)%buoy_lat, buoy_coll(kbuoy)%buoy_lon, &
          trim(dtg1), buoy_speed, buoy_coll(kbuoy)%buoy_dir, &
          buoy_coll(kbuoy)%scat_lat, buoy_coll(kbuoy)%scat_lon, trim(dtg2), &
          buoy_coll(kbuoy)%scat_speed, buoy_coll(kbuoy)%scat_dir, &
          buoy_coll(kbuoy)%model_speed, buoy_coll(kbuoy)%model_dir, &
          buoy_coll(kbuoy)%wvc, buoy_coll(kbuoy)%alfa,  buoy_coll(kbuoy)%mle
      end if
    endif
  enddo
  !  #]

  ! close input and output files
  if (use_bufr_scat_input) then
     call close_scat_bufr_file()
  end if
  if (use_ascii_scat_input) then
     call close_scat_ascii_file()
  end if

  close(unit_ascii_output)

  if (allocated(buoy_coll)) deallocate(buoy_coll)

  ! format statements
  101 format (I4.4,2I2.2,1X,2I2.2)
  102 format (I5,F7.2,F8.2,2X,A,2F6.1,F9.2,F8.2,2X,A,4F6.1,I4,2F8.2)
  103 format (I5,F7.2,F8.2,2X,A,2F6.1,F9.2,F8.2,2X,A,4F6.1,I4,2F8.2,2f6.1)

contains
  !---------------------------------
  subroutine open_scat_bufr_file(filename_bufr)
    !  #[ open the file and init loop counters
    character(len=256), intent(in) :: filename_bufr

    ! open BUFR input file and obtain number of messages present
    call set_BUFR_fileattributes(filename_bufr, bufr_file_attributes, &
                                 error_flag)
    call open_BUFR_file(bufr_file_attributes, error_flag)
    nr_of_messages = get_BUFR_nr_of_messages(bufr_file_attributes)
    
    ! read first BUFR message and retrieve some information from it
    call get_BUFR_message(bufr_file_attributes, 1, tbd, error_flag)
    
    ! choose which BUFR template to use when interpreting the data
    ! (i.e. determine type of message depending on BUFR descriptor)
    ! this assumes all bufr messages in the file use the same template.
    call set_BUFR_template(tbd)
    
    kmsg = 0
    ksub = 0

  end subroutine open_scat_bufr_file
    !  #]
  subroutine get_next_bufr_scat_result(res, wvc, lat, lon, mspd, mdir, qcflag, &
       alfa, aspd, adir, spd, dir, mle, year, month, day, hour, minute, &
       research_option_set, end_of_file_reached)
    !  #[ get next scat wvc result from file
    integer, intent(out) :: res, wvc
    real,    intent(out) :: lat, lon, mspd, mdir
    logical, intent(out) :: qcflag
    real,    intent(out) :: alfa, aspd, adir, spd, dir, mle
    integer, intent(out) :: year, month, day, hour, minute
    logical, intent(out) :: research_option_set, end_of_file_reached

    if ( (kmsg .eq. 0) .or. &
         ((kmsg .gt. 0) .and. (ksub .eq. tbd%nsubsets)) ) then
       ! start
       kmsg = kmsg+1

       if (kmsg .gt. nr_of_messages) then
          end_of_file_reached = .true.
          return
       end if

       call get_BUFR_message(bufr_file_attributes, kmsg, tbd, error_flag)
       
       ! check if BUFR file was processed with -research option. If so,
       ! it will contain the 2DVAR analysis wind speed and direction 
       ! instead of some radar parameters.
       
       if (get_sw_id_l1(tbd) == 0) then
          research_option_set = .true.
          if (kmsg .eq. 1) print '(A)','Data contains 2DVAR analysis wind'
       end if

    end if

    if (ksub .lt. tbd%nsubsets) then
       ksub = ksub+1
    else
       ksub = 1
    end if

    ! step to the next subset
    call goto_subset(ksub)

    res  = get_res(tbd)
    lat  = get_lat(tbd)
    lon  = get_lon(tbd)
    mspd = get_mspd(tbd)
    mdir = get_mdir(tbd)
    qcflag  = get_qcflag(tbd)
    alfa = get_trk(tbd)
    wvc  = get_cln(tbd)
    
    aspd = -999.9
    adir=  -999.9
    if (research_option_set) then
       aspd = get_aspd(tbd)
       adir = get_adir(tbd)
    end if

    spd = get_spd(tbd)
    dir = get_dir(tbd)

    ! NOTE: this is the MLE of the first rank solution
    ! use closest=.false. for the MLE of the selected solution
    mle = get_bsdist(tbd,closest=.true.)

    year   = get_year(tbd)
    month  = get_month(tbd)
    day    = get_day(tbd)
    hour   = get_hour(tbd)
    minute = get_minute(tbd)

  end subroutine get_next_bufr_scat_result
    !  #]
  subroutine close_scat_bufr_file
    !  #[ close the BUFR file

    call close_BUFR_file(bufr_file_attributes, error_flag)

  end subroutine close_scat_bufr_file
    !  #]
  !---------------------------------
  subroutine open_scat_ascii_file(filename_ascii, instrument)
    !  #[ open the file
    character(len=256), intent(in) :: filename_ascii
    integer, intent(in) :: instrument
    
    ! local variables
    integer :: status
    character(len=256) :: tmpstr

    current_instrument = instrument

    open(unit_scat_ascii_input, file=trim(filename_ascii), status='old')

    ! read res from file header
    read(unit_scat_ascii_input, fmt='(A)', iostat=status, end=999) tmpstr
    if (status .ne. 0) then
       print *,'ERROR: could not read first header line'
    end if

    read(unit_scat_ascii_input, fmt='(A)', iostat=status, end=999) tmpstr
    if (status .ne. 0) then
       print *,'ERROR: could not read second header line'
    end if

    !print *,'tmpstr = '//trim(tmpstr)
    !print *,'tmpstr(8:) = ['//trim(tmpstr(8:))//']'
    !tmpstr = trim(tmpstr(8:)
    read(tmpstr(8:),*) stored_res

    ! expected format:
    !
    !#lat, lon, year, month, day, hour, minute, w_spd, w_dir, mle, qcflag, m_spd, m_dir, wvc
    !# res = 25000
    !4.23 310.84 2012  1  1  0 57   6.88 235.90   0.00       0 999.00 999.00  22
    !4.28 311.06 2012  1  1  0 57   7.22 235.80  -0.20       0 999.00 999.00  23
    !4.32 311.28 2012  1  1  0 57   7.61 234.70  -0.20       0 999.00 999.00  24

    return

999 print *,'ERROR: unexpected end of file...'
    stop

  end subroutine open_scat_ascii_file
    !  #]
  subroutine get_next_ascii_scat_result(res, wvc, lat, lon, mspd, mdir, qcflag,&
       alfa, aspd, adir, spd, dir, mle, year, month, day, hour, minute, &
       research_option_set, end_of_file_reached)
    !  #[ get next scat wvc result from file
    integer, intent(out) :: res, wvc
    real,    intent(out) :: lat, lon, mspd, mdir
    logical, intent(out) :: qcflag
    real,    intent(out) :: alfa, aspd, adir, spd, dir, mle
    integer, intent(out) :: year, month, day, hour, minute
    logical, intent(out) :: research_option_set, end_of_file_reached

    ! local variables
    integer :: status, qcflag_int
    logical :: data_found
    character(len=256) :: tmpstr

    res    = stored_res
    wvc    = -1
    lat    = 999.
    lon    = 999.
    mspd   = 999.
    mdir   = 999.
    qcflag_int = -1
    qcflag = .false.
    alfa   = 999.
    aspd   = 999.
    adir   = 999.
    spd    = 999.
    dir    = 999.
    mle    = 999.
    year   = -1
    month  = -1
    day    = -1
    hour   = -1
    minute = -1
    research_option_set = .false.

    error_flag = no_error

    ! expected ascii data format:
    !
    !#lat, lon, year, month, day, hour, minute w_spd, w_dir, mle, qcflag, m_spd, m_dir, wvc
    !# res = 25000
    !  6.31 306.22 2012  6  2  1 33   6.26 248.80  -0.10   32768   7.80 252.30  42
    !  6.43 305.73 2012  6  2  1 33   5.80 249.50   0.00   32768   7.90 252.40  40
    !  6.48 305.95 2012  6  2  1 33   6.08 249.80   0.00   32768   7.88 252.20  41

    data_found = .false.
    do while (.not. data_found)
       read(unit_scat_ascii_input, fmt='(A)', iostat=status, end=999) tmpstr
       if ( (status .eq. 0) .and. &
            (tmpstr(1:1) .ne. '#')           ) then
          read(tmpstr,*, iostat=status, end=999) lat, lon, year, month, day, hour, minute, &
               spd, dir, mle, qcflag_int, mspd, mdir, wvc
          if (status .ne. 0) then
             print *,' '
             print *,'File '//trim(filename_buoy)//&
                  ' seems to contain invalid buoy data near line', nbuoys
             print *,' '
             print *,'Skipping this dataline:  ['//trim(tmpstr)//']'
             print *,' '
          endif
          data_found = .true.
       endif
    enddo
    
    qcflag = qc_int_to_flag(qcflag_int)

    end_of_file_reached = .false.
    if (.not. data_found) end_of_file_reached = .true.

    !print *, 'DEBUG:', lat, lon, year, month, day, hour, minute, &
    !           spd, dir, mle, qcflag_int, mspd, mdir, wvc

    return

999 end_of_file_reached = .true.
    return

  end subroutine get_next_ascii_scat_result
    !  #]
  function qc_int_to_flag(qcflag_int) result(qcflag)
    !  #[ convert qcflag integer into a boolean
    integer, intent(in) :: qcflag_int
    logical :: qcflag ! result

    ! local variables                                                           
    integer :: bit1_to_inspect, bit2_to_inspect

    bit1_to_inspect = -1
    bit2_to_inspect = -1
    if (current_instrument .eq. instr_is_seawinds) then
       bit1_to_inspect = 10 ! KNMI QC flag                                      
       bit2_to_inspect = 2  ! NOAA rain flag                                    
    else ! valid for GenSeawinds and ascat                                      
       bit1_to_inspect = 17 ! KNMI QC flag                                      
    end if

    if (missing_int(qcflag_int)) then
       qcflag = .true.
       return
    else
       ! do not use data for which QC flag is set                               
       if (bit1_to_inspect .ne. -1) then
          qcflag = btest(qcflag_int, bit1_to_inspect)
          if (qcflag) return
       end if
       if (bit2_to_inspect .ne. -1) then
          qcflag = btest(qcflag_int, bit2_to_inspect)
          if (qcflag) return
       end if
    endif

  end function qc_int_to_flag
    !  #]
  subroutine close_scat_ascii_file
    !  #[ close the file
     close(unit_scat_ascii_input)

  end subroutine close_scat_ascii_file
    !  #]
  !---------------------------------
end program BuoyCollocate

